# -*- coding: utf-8 -*-
"""
Created on Thu May 11 15:23:51 2023

@author: lenovo
"""

import numpy as np
import pandas as pd
from osgeo import gdal


# import os
# import glob
# import shutil
# from tqdm import tqdm # 运行计时

################################################
def read_img(filename):
    dataset = gdal.Open(filename)  # 打开文件
    if dataset == None:
        print(filename + "文件无法打开")
        return
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    # im_bands = dataset.RasterCount #波段数
    # im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵
    # im_proj = dataset.GetProjection()  # 地图投影信息
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height).astype(np.float32)  # 将数据写成数组，对应栅格矩阵
    # im_lon=[im_geotrans[0]+i*im_geotrans[1] for i in range(im_width)]
    # im_lat=[im_geotrans[3]+i*im_geotrans[5] for i in range(im_height)]

    del dataset  # 关闭对象，文件dataset   
    return im_data


# =========================
#
town = ['商州区', '洛南县', '丹凤县', '商南县', '山阳县', '镇安县', '柞水县', '商洛市']

# excel表格的行标签
# vartype = ['水源涵养', '土壤保持', '洪水调蓄',
#            '空气净化', '水质净化', '固碳价值',
#            '释氧价值', '气候调节', '负氧离子',
#            '储量碳价值', '水资源']

vartype = [
    '水源涵养', '土壤保持', '洪水调蓄',
    '空气净化', '水质净化', '固碳价值',
    '释氧价值', '气候调节', '负氧离子',
]

# tifpath=r'H:\company\公司项目\陕西\商洛\data\masked\Alberts'
# df=pd.read_excel(excelpath)


# 保存的结果路径
excelpath = r'H:\company\公司项目\陕西\商洛\data\统计结果'
fntt = excelpath + '\\总量.xlsx'
fnav = excelpath + '\\平均.xlsx'

# 要保存到excel里的数据
dft = pd.DataFrame(columns=['生态价值(亿元)'], data=vartype)
dfa = pd.DataFrame(columns=['生态价值(元/m2)'], data=vartype)
# ci遍历的城镇名字town
for ci in town:
    # 统计的值
    tt = []
    # 平均的值
    aa = []
    for tp in vartype:
        # tif文件的地址
        # tifpath=r'H:\company\公司项目\陕西\商洛\data\masked\Alberts'
        tifpath = r'F:\公司项目\陕西\商洛\data\masked\Alberts'
        # fin读取的tif文件路径 示例：F:\公司项目\陕西\商洛\data\masked\Alberts\水源涵养_商洛市_Alberts.tif
        fin = tifpath + '\\' + tp + '_' + ci + '_Alberts.tif'
        print('fin', fin)
        im_data = read_img(fin)
        im_data[im_data == -9999] = np.nan
        tot1 = np.nansum(im_data)
        tot2 = np.nansum(tot1)
        tot = tot2 * 30 * 30
        tot = tot / 10000 / 10000
        ave = np.nanmean(np.nanmean(im_data))
        tt.append(tot)
        aa.append(ave)
    dft[ci] = tt
    dfa[ci] = aa

dft.to_excel(fntt)
dfa.to_excel(fnav)

# intif = r'H:\company\公司项目\河南\新野县\插值结果\TiffFile\Alberts\nep_2011_Alberts.tif'

# im_data=read_img(intif)

# im_data[im_data==-99]=np.nan
# im_data[~np.isnan(im_data)]=1
# area=np.nansum(im_data)*100/1e6
# print(area)
